Data wrangel

data(UCBadmit, package = "rethinking")
d <- UCBadmit %>% 
  mutate(gid = ifelse(applicant.gender == "male", "1", "2"))
rm(UCBadmit)

d
##    dept applicant.gender admit reject applications gid
## 1     A             male   512    313          825   1
## 2     A           female    89     19          108   2
## 3     B             male   353    207          560   1
## 4     B           female    17      8           25   2
## 5     C             male   120    205          325   1
## 6     C           female   202    391          593   2
## 7     D             male   138    279          417   1
## 8     D           female   131    244          375   2
## 9     E             male    53    138          191   1
## 10    E           female    94    299          393   2
## 11    F             male    22    351          373   1
## 12    F           female    24    317          341   2

# models

b11.7 <-
  brm(data = d, 
      family = binomial,
      admit | trials(applications) ~ 0 + gid,
      prior(normal(0, 1.5), class = b),
      iter = 2000, warmup = 1000, cores = 4, chains = 4,
      seed = 42) 

b11.8 <-
  brm(data = d, 
      family = binomial,
      bf(admit | trials(applications) ~ a + d,
         a ~ 0 + gid, 
         d ~ 0 + dept,
         nl = TRUE),
      prior = c(prior(normal(0, 1.5), nlpar = a),
                prior(normal(0, 1.5), nlpar = d)),
      iter = 4000, warmup = 1000, cores = 4, chains = 4,
      seed = 42) 

# Creating custom distribution for m12.1

beta_binomial2 <- custom_family(
  "beta_binomial2", dpars = c("mu", "phi"),
  links = c("logit", "log"), lb = c(NA, 2),
  type = "int", vars = "vint1[n]"
)

stan_funs <- "
  real beta_binomial2_lpmf(int y, real mu, real phi, int T) {
    return beta_binomial_lpmf(y | T, mu * phi, (1 - mu) * phi);
  }
  int beta_binomial2_rng(real mu, real phi, int T) {
    return beta_binomial_rng(T, mu * phi, (1 - mu) * phi);
  }"

stanvars <- stanvar(scode = stan_funs, block = "functions")

b12.1 <-
  brm(data = d, 
      family = beta_binomial2,  # the custom likelihood
      admit | vint(applications) ~ 0 + gid,
      prior = c(prior(normal(0, 1.5), class = b),
                prior(exponential(1), class = phi)),
      iter = 2000, warmup = 1000, cores = 4, chains = 4,
      stanvars = stanvars,  # note our `stanvars`
      seed = 42)


b11.7_dept <- brm(data = d, 
      family = binomial,
      admit | trials(applications) ~ 0 + gid + (0 + gid | dept),
      prior(normal(0, 1.5), class = b),
      prior(exponential(1), class = sd, group = dept), 
      prior(lkj(2), class = cor, group = dept),
        
      iter = 2000, warmup = 1000, cores = 4, chains = 4,
      seed = 42) 
models <- list(b11.7, b11.8, b12.1, b11.7_dept)
pp_models <- list(b11.7, b11.8, b11.7_dept)

purrr::map(pp_models, pp_check)
## Using 10 posterior samples for ppc type 'dens_overlay' by default.
## Using 10 posterior samples for ppc type 'dens_overlay' by default.
## Using 10 posterior samples for ppc type 'dens_overlay' by default.
## [[1]]

## 
## [[2]]

## 
## [[3]]

purrr::map(models, plot)

## [[1]]
## [[1]][[1]]

## 
## 
## [[2]]
## [[2]][[1]]

## 
## [[2]][[2]]

## 
## 
## [[3]]
## [[3]][[1]]

## 
## 
## [[4]]
## [[4]][[1]]

purrr::map(models, base::summary)
## [[1]]
##  Family: binomial 
##   Links: mu = logit 
## Formula: admit | trials(applications) ~ 0 + gid 
##    Data: d (Number of observations: 12) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## gid1    -0.22      0.04    -0.30    -0.15 1.00     3753     2680
## gid2    -0.83      0.05    -0.93    -0.73 1.00     2663     2449
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## 
## [[2]]
##  Family: binomial 
##   Links: mu = logit 
## Formula: admit | trials(applications) ~ a + d 
##          a ~ 0 + gid
##          d ~ 0 + dept
##    Data: d (Number of observations: 12) 
## Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
##          total post-warmup samples = 12000
## 
## Population-Level Effects: 
##         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## a_gid1     -0.53      0.53    -1.61     0.54 1.00      982     1196
## a_gid2     -0.43      0.54    -1.51     0.65 1.00      986     1238
## d_deptA     1.11      0.54     0.02     2.20 1.00      990     1232
## d_deptB     1.06      0.54    -0.02     2.15 1.00      995     1226
## d_deptC    -0.15      0.54    -1.23     0.94 1.00      987     1266
## d_deptD    -0.19      0.54    -1.26     0.91 1.00      991     1210
## d_deptE    -0.63      0.54    -1.71     0.46 1.00     1004     1267
## d_deptF    -2.19      0.55    -3.29    -1.08 1.00     1023     1336
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## 
## [[3]]
##  Family: beta_binomial2 
##   Links: mu = logit; phi = identity 
## Formula: admit | vint(applications) ~ 0 + gid 
##    Data: d (Number of observations: 12) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## gid1    -0.44      0.42    -1.25     0.37 1.00     3009     2283
## gid2    -0.32      0.44    -1.18     0.55 1.00     3301     2277
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi     3.01      0.78     2.04     4.94 1.00     1909     1153
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## 
## [[4]]
##  Family: binomial 
##   Links: mu = logit 
## Formula: admit | trials(applications) ~ 0 + gid + (0 + gid | dept) 
##          autocor ~ tructure(list(), class = "formula", .Environment = <environment>)
##    Data: d (Number of observations: 12) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~dept (Number of levels: 6) 
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(gid1)           1.28      0.41     0.72     2.30 1.00     1279     1652
## sd(gid2)           1.55      0.49     0.88     2.77 1.00     1332     1735
## cor(gid1,gid2)     0.85      0.18     0.33     0.99 1.00     1102     2196
## 
## Population-Level Effects: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## gid1    -0.52      0.47    -1.46     0.49 1.00     1140     1761
## gid2    -0.32      0.57    -1.44     0.82 1.00     1404     2112
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
purrr::map(models, mcmc_plot)
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

expose_functions(b12.1, vectorize = TRUE)
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Users/shawn/Library/R/4.0/library/Rcpp/include/"  -I"/Users/shawn/Library/R/4.0/library/RcppEigen/include/"  -I"/Users/shawn/Library/R/4.0/library/RcppEigen/include/unsupported"  -I"/Users/shawn/Library/R/4.0/library/BH/include" -I"/Users/shawn/Library/R/4.0/library/StanHeaders/include/src/"  -I"/Users/shawn/Library/R/4.0/library/StanHeaders/include/"  -I"/Users/shawn/Library/R/4.0/library/RcppParallel/include/"  -I"/Users/shawn/Library/R/4.0/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Users/shawn/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Users/shawn/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Users/shawn/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Users/shawn/Library/R/4.0/library/RcppEigen/include/Eigen/Core:88:
## /Users/shawn/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Users/shawn/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Users/shawn/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Users/shawn/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
## /Users/shawn/Library/R/4.0/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
# "Unfortunately, using a custom distribution means that we can’t easily let STAN compute the ELPD for us, so we have to do it manually in R"
log_lik_beta_binomial2 <- function(i, draws) {
  mu <- draws$dpars$mu[, i]
  phi <- draws$dpars$phi
  trials <- draws$data$vint1[i]
  y <- draws$data$Y[i]
  beta_binomial2_lpmf(y, mu, phi, trials)
}

# loo compair of all models

loo(b11.7, b11.8, b12.1, b11.7_dept)
## Warning: Found 6 observations with a pareto_k > 0.7 in model 'b11.7'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## Warning: Found 5 observations with a pareto_k > 0.7 in model 'b11.8'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## Warning: Found 11 observations with a pareto_k > 0.7 in model 'b11.7_dept'. It
## is recommended to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
## Warning: Not all models have the same y variable. ('yhash' attributes do not
## match)
## Output of model 'b11.7':
## 
## Computed from 4000 by 12 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo   -480.1 157.4
## p_loo       100.4  34.0
## looic       960.2 314.9
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     3     25.0%   613       
##  (0.5, 0.7]   (ok)       3     25.0%   172       
##    (0.7, 1]   (bad)      1      8.3%   35        
##    (1, Inf)   (very bad) 5     41.7%   1         
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'b11.8':
## 
## Computed from 12000 by 12 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo    -57.3  8.6
## p_loo        12.4  3.6
## looic       114.6 17.1
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     3     25.0%   2808      
##  (0.5, 0.7]   (ok)       4     33.3%   267       
##    (0.7, 1]   (bad)      3     25.0%   146       
##    (1, Inf)   (very bad) 2     16.7%   18        
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'b12.1':
## 
## Computed from 4000 by 12 log-likelihood matrix
## 
##          Estimate  SE
## elpd_loo    -68.9 3.0
## p_loo         2.0 0.6
## looic       137.8 6.0
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     11    91.7%   1112      
##  (0.5, 0.7]   (ok)        1     8.3%   990       
##    (0.7, 1]   (bad)       0     0.0%   <NA>      
##    (1, Inf)   (very bad)  0     0.0%   <NA>      
## 
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'b11.7_dept':
## 
## Computed from 4000 by 12 log-likelihood matrix
## 
##          Estimate  SE
## elpd_loo    -49.7 2.2
## p_loo        11.1 1.2
## looic        99.5 4.4
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     0      0.0%   <NA>      
##  (0.5, 0.7]   (ok)       1      8.3%   299       
##    (0.7, 1]   (bad)      6     50.0%   58        
##    (1, Inf)   (very bad) 5     41.7%   11        
## See help('pareto-k-diagnostic') for details.
## 
## Model comparisons:
##            elpd_diff se_diff
## b11.7_dept    0.0       0.0 
## b11.8        -7.6       7.8 
## b12.1       -19.2       1.4 
## b11.7      -430.4     156.9

Using the basic loo output we see that m11.7_dept has the lowest lepd. Interesting here is that its version b11.7 was the worst out of the 4 models, and by very large margin.

# previous method of waic and loo compare no longer seem to work due to the custom distribution of m12.1. The commented-out code is included below 
waic_b11.7 <- add_criterion(b11.7, criterion = "waic")
waic_b11.8 <- add_criterion(b11.8, criterion = "waic")
waic_b12.1 <- add_criterion(b12.1, criterion = "waic")
waic_b11.7_dept <- add_criterion(b11.7_dept, criterion = "waic")

comparison_waic <- loo_compare(waic_b11.7, waic_b11.8, waic_b12.1, waic_b11.7_dept, criterion = "waic")
print(comparison_waic, simplify = FALSE)

loo_b11.7 <- add_criterion(b11.7, criterion = "loo")
loo_b11.8 <- add_criterion(b11.8, criterion = "loo")

expose_functions(b12.1, vectorize = TRUE)

# "Unfortunately, using a custom distribution means that we can’t easily let STAN compute the ELPD for us, so we have to do it manually in R"
log_lik_beta_binomial2 <- function(i, draws) {
  mu <- draws$dpars$mu[, i]
  phi <- draws$dpars$phi
  trials <- draws$data$vint1[i]
  y <- draws$data$Y[i]
  beta_binomial2_lpmf(y, mu, phi, trials)
}

loo_b12.1 <- loo(b12.1)

loo_b12.1 <- add_criterion(b12.1, criterion = "loo")
loo_b11.7_dept <- add_criterion(b11.7_dept, criterion = "loo")


comparison_loo <- loo_compare(loo_b11.7, loo_b12.1, b12.1, loo_b11.7_dept, criterion = "loo")
print(comparison_loo, simplify = FALSE)